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ABSTRACT 

Core radii of globular clusters in the Large and Small Magellanic Clouds show an increasing trend with age. 
We propose that this trend is a dynamical effect resulting from the accumulation of massive stars and stellar- 
mass black holes at the cluster centers. The black holes are remnants of stars with initial masses exceeding 
~ 2O-25M0; as their orbits decay by dynamical friction, they heat the stellar background and create a core. 
Using analytical estimates and A^-body experiments, we show that the sizes of the cores so produced and their 
growth rates are consistent with what is observed. We propose that this mechanism is responsible for the 
formation of cores in all globular clusters and possibly in other systems as well. 
Subject headings: black hole physics — gravitation — gravitational waves — galaxies: nuclei 



1. INTRODUCTION 

An enduring problem is the origin of cores, regions near the 
center of a stellar or dark matter system where the density is 
nearly constant. Resolved cores clearly exist in some stellar 
systems, e. g. globular clusters ( Harris 1996). In other sys- 
tems, such as early-type galaxies, cores were long believed to 
be generic but we re later shown to be artifacts of the seeing 
(|§chweizer Il979h . Ne vertheless a few el liptical galaxies do 
exhibit bona-fide cores JLauer et al. F20Q2) while many others 
show a central density that rises only very slowly toward the 
center ( Merritt & Fridman 1995). Density profiles of struc- 
tures that form from gravitational clustering of density per- 
turbations in an expanding universe are believed to lack cores 
(|Power et al."2003), although there is evidence for dark mat- 
ter cores in the rotation curves of some late-type galaxies (e.g. 
tlimenez. Verde & OM l l2Q03h ). 

The existence of a core is usually deemed to require a spe- 
cial explanation. For instance, galaxy cores may form when 
binary black holes eject stars via the gravitational slingshot 
(pibisuzaki, Makino & Okumura 1991). 

A useful sample for testing theories of core formation is 
the ensemble of globular clusters (GCs) around the Large and 
Small Magellanic Clouds (LMC/SMC). These clusters have 
masses similar to those of Galactic GCs, but many are much 
younger, with ag es that range from 10^ - lO'" yr Furthermore 
ground-based ( Elson, Freeman & Lauer 1989; Elson 199 J, 
Il992) and HST (Mackev & Gilmore 2003 a b) observations 
reveal a clear trend of core radius with age: while young 
clusters (r ^ 10^ yr) have core radii consistent with zero, 
clusters older than ~ 10^ yr exhibit the full range of core 
sizes seen in Galactic GCs, pc < < 10 pc (Figure 
The maximum core radius observed in the LMC/SMC GCs 
is an increasing function of age and is given roughly by 
rc « 2.25pclogjoTyr- 14.5. Attempts to explain the core ra- 
dius evolution in terms of stellar mass loss (Elson 1991), ^ 
primordial population of binary stars, or time-varying tidal 
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fields (^Wilk inson et a l."2003') have met with limited success. 
The difficulty is to find a mechanism that can produce sub- 
stantial changes in the central structure of a GC on time scales 
as short as a few hundred Myr, while leaving the large-scale 
structure of the cluster intact. 

In this paper, we describe a new mechanism for the forma- 
tion of GC cores and their evolution with time. Massive stars 
and their black hole remnants sink to the center of a GC due 
to dynamical friction against the less massive stars. The en- 
ergy transferred to the stars during this process, and during 
the three- and higher-A^ encounters between the black holes 
that follow, has the effect of displacing the stars and creating 
a core. The rate of core growth implied by this model is con- 
sistent with the observed dependence of core size on age in 
the LMC/SMC clusters. 
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Fig. 1. — Core radius versus age for LMC and SMC GCs from the 
samples of Mackey & Gilmore (2003a,b). Lines show core radius evolu- 
tion from the W-body simulations with initial cusp slope 7 = 1 and three 
different scaUngs to physical units; see text for details. 



2. CORE FORMATION TIMESCALES 

Consider a gravitationally bound stellar system in which 
most of the mass is in the form of stars of mass m, but which 
also contains a subpopulation of more massive objects with 
masses niBH- The orbits of the more massive objects de- 
cay due to dynamical friction. Assume that the stellar den- 
sity profile is initially a power-law in radius, p{r) = K(r/a)~^, 
K = (3-7)M/47rfl-' with M the total stellar mass and a the 
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density scale length; the expression for K assumes that the 
density follows a Dehnen (1993) law outside of the stellar 
cusp, i.e. p ~ r""* at large r. The effective radius (the radius 
containing 1 /2 of the mass in projection) is related to a via 
7;,/fl« (1.8, 1.3, 1.0) for 7 = (1,1.5,2). 

Due to the high central concentration of the mass, the orbits 
of the massive particles will rapidly circularize as they receive 
nearly-impulsive velocity changes near pericenter. Once cir- 
cular, orbits shrink at a rate that can be computed by equating 
the torque from dynamical friction with the rate of change of 
orbit al angular mo mentum. We adopt the usual approxima- 
tion jSr)it7eHI1987h in which the frictional force is produced 
by stars with velocities less than the orbital velocity of the 
massive object. The rate of change of the orbital radius, as- 
suming a fixed and isotropic stellar background, is then 



dt 4-7 \ a M \a 
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where /3 = (6-7)/2(2-7) and In A is the Coulomb logarithm, 
roughly equal to 6.6 (Spinnato, Fellhauer & Porteaies Zwart 
Il003). For7 = (1.0,1.5,2.0),F = (0.193, 0.302, 0.427). Equa- 
tion (1) implies that the massive object comes to rest at the 
center of the stellar system in a time 
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with r,- the initial orbital radius; the leading coefficient de- 
pends weakly on 7. Equation (|2j can be written 



Af w3x Wyra^M. 
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with ajo the density scale length in units of 10 pc (e.g. Figure 
1 of Ivan den Bergh (1991)), M5 = M/IQ^Mq, and mBH.m = 
mBH/iOMQ, the approximate masses of black hole r emnants 
of stars with initial masses exceeding ~ 2 - 25 M(r, (iMaederi 
[1992; Portegies Zwart, Verbunt & Ergml ll997h . This time 
is of the same order as the time (^ 10^ yr) over which 
core expansion is observed to take place ( iMackev & Gilmora 
(|2003ab). Fig. 1). 

To estimate the effect of the massive remnants on the stel- 
lar density profile, consider the evolution of an ensemble of 
massive particles in a stellar system with initial density pro- 
file p ^ r~^. The energy released as one particle spirals in 
from radius r,- to rj is 2mBH<^^in{rj/rf), with <t the ID stellar 
velocity dispersion. Decay will halt when the massive par- 
ticles form a self-gravitating system of radius ^ GMbh/o''^ 
with Mbh = X]'"BH- Equating the energy released during in- 
fall with the energy of the stellar matter initially within r^, the 
"core radius," gives 



2GMbh f n<j^ 



\GM, 



BH 



(4) 



Most of the massive particles that deposit their energy 
within will come from radii r,- « a few x r^, im- 
plying w several x GMbh/o'^ and a displaced stel- 
lar mass of ~ several x Mrh- If Mbh ~ 10"^M 
iPortegies Zwart & McMillanll2000h . then Vc/a w several x 



2Mbh/M and the core radius is roughly 10% of the effective 
radius. 

Evolution continues as the massive particles form binaries 
and begin to engage in three-body interactions with other mas- 
sive particles. These superelastic encounters will eventually 
eject most or all of the massive particles from the cluster 
Assume that this ejection occurs via the cumulative effect 
of many encounters, such that almost all of the binding en- 
ergy so released can find its way into the stellar system as 
the particle spirals back into the core. The energy released 
by a single binary in shrinking to a separation such that its 
orbital velocity equals the escape velocity from the core is 
~ mBHO'^in(4MBH/M). If all of the massive particles find 
themselves in such binaries before their final ejection and if 
most of their energy is deposited near the center of the stellar 
system, the additional core mass will be 



, M \ 
Mbh J 



(5) 



e.g. ~ 5Mbh for M /Mbh = 100, similar to the mass displaced 
by the initial infall. The additional mass displacement takes 
place over a much longer time scale however and additional 
processes (e.g. core collapse) may compete with it. 
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Fig. 2. — Evolution of the mass deficit, from the W-body experiments 
with Nbh = 10. Data were smoothed with cubic spHnes. 

3. W-BODY SIMULATIONS 

We used A^-body simulations to test the core formation 
mechanism described above. Initial conditions were designed 
to represent GCs in which 1 % of the total mass is initially in 
the form of massive objects, either stars or their black hole 
remnants. Integrations were carried out using NBODY6-H-, 
a high-precision, parallel, fourth-order direct force integrator 
which implements coordinate regularization for close encoun- 
ters (Spurzem 1999). Particles had one of two masses, repre- 
senting either black holes {mbh) or stars {m). The number of 
particles representing black holes was Nbh = (4, 10,20) and 
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the ratio of ntBH to m ranged from 10 to 25. Most of the A^- 
body experiments used = lO'* particles. We concentrate here 
on the results obtained withNBH= 10 and mBH/in= 10; results 
obtained with other values of Nbh were consistent. All parti- 
cles were initially distributed according to Dehnen's (1993) 
density law, p(r) = ((3 - -f)M/4na^)/C /(I + 0^'^ ,C = r/a. 
The logarithmic slope of the central density cusp is specified 
via the parameter 7. Initial velocities were generated assum- 
ing isotropy; positions and velocities of the massive particles 
were distributed in the same way as the stars, so that the initial 
conditions represented a cluster in which the massive objects 
had not yet begun to segregate spatially with respect to the 
lighter stars. 

The 1% mass fraction in massive particles was based on 
a Scalo (1986) mass distribution with lower and upper mass 
limits of 0.1 Mq and 100 Mq respectively. With such a mass 
function, about 0.071% of the stars are more massive than 
2OM0 and 0.045% are more massive than 25M0. A star 
cluster containing A^^ stars thus produces ~ 6 x 10~''A^* black 
holes. Known Galactic black holes have masses mbh be- 
tween 6M0 and I8M0 (Timmes, Wooslev, & Weaver 199dl) . 
Adopting an average black hole mass of IOM0 then results in 
a total black hole mass of ~6x lO-^M. 

The decision to use just two mass groups - clearly an 
idealization of the true situation - was made for two rea- 
sons. First, the interpretation of the A^-body results is greatly 
simplified in such a model. Second, it is not clear what a 
better choice for the initial spectrum of masses would be. 
The distribution of blac k hole masses produced by stars wit h 
m > 25Mq is uncertain (Timmes. Wooslev. & Weaveil 19961) : 
some Gala ctic black holes may have masses as low as ^ 3M0 
JWhite & va n Paradiis 1996), close to the maximum probable 
masses of neutron stars. However this may be a selection ef- 
fect: low-mass binary systems tend to be selected due to their 
longer lifetimes. The remnant mass range between ^ IM0 
and ^ 3M0 is occupied by neutron stars, but there is evi- 
dence that neutron sta rs receive larger kicks at birth than black 
holes (iLvne & Lorimer..l994) and may be ejected. (Portegies 
Zwart & McMillan (2000) estimate that ~ 10% of black holes 
are ejected from GCs by formation kicks.) In summary, the 
initial spectrum of masses in a GC shortly after its forma- 
tion is poorly known. Future studies will attempt to include 
a realistic treatment of stellar physics, primordial binaries, 
star formation, and other processes that affect the initial mass 
spectrum and the initial spatial distribution of different mass 
groups. 

If the core radius is defined as the radius at which the pro- 
jected density falls to 1/2 of its central value, Dehnen models 
with 7 > 1 have rc = 0. Any core that appears in these models 
must therefore be a result of dynamical evolution. 

Henceforth we adopt units in which G = a = M = 1 . The 
corresponding unit of time is 
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where Mg is the cluster mass in units of 1O^M0 and a 10 is 
the cluster scale length in units of 10 pc. The effective radius 
Re, defined as the radius containing 1 /2 of the light particles 
seen in projection, is « (1.8, 1.3, 1.0) in model units (a = 
1) for 7 = (1, 1.5,2). The time scaling of equation is not 
correct for processes whose rates depend on the masses of 
individual stars or black holes, since our models have fewer 
stars than real GCs. The most important of these processes 



for our purposes are black hole-star interactions, which are 
responsible for the orbital decay of the black holes and the 
growth of the core. This decay occurs in our simulations at a 
rate that is ^ N,/Nbh times faster than implied by the scaling 
of equation (|6j, with A^, the true number of black holes in 
a GC. Assuming a Scalo initial mass function as above, this 
factor may be written 6.0(A^*/10^)/(A^bh/10) with A^* the 
true number of stars in a GC. 

As discussed above, we expect the stellar mass displaced 
by the massive particles to scale roughly with Mbh- One 
way to illustrate this is via the mass deficit, defined as in 
Milosavlievic et al. (2002): it is the mass difference between 
the initial stellar density /9(r,0) and the density at time t, in- 
tegrated from the origin out to the radius at which p(r, t) first 
exceeds p(r, 0). The mass deficit is a measure of the core mass. 
Figure 2 shows M^siit)/ Mbh for the A^-body experiments. The 
density center was computed via the Casertano-Hut (1985) al- 
gorithm. The black holes displace a mass in stars of order 
2-8 times their own mass; the larger values correspond to the 
larger values of 7 although there is considerable scatter from 
experiment to experiment for a given 7. The results for 7 = 2 
are consistent with the analytic arguments presented above, 
which implied a core mass of a few times Mbh after a time in 
model units of ^ Q.2M/mBH ~ 20 (cf. equation|2ji followed 
by a slower displacement of a similar mass as the black hole 
particles engage in three-body interactions (cf. equation|5}. 

Some of the A^-body simulations show a decrease in the 
core radius after t « lO"* (Figure 2, 3). By this time, the ma- 
jority of the black holes have been ejected. Our simulated 
clusters are then effectively reduced to equal-mass systems, 
which take about 15 half mass relaxation times to experience 
core collapse (Spitzer 1987). The two-body relaxation time is 
roughly Tr w Q.ITqN / \viN w 200 in our A^-body models, with 
To the crossing time. It is therefore not surprising that, once 
the black holes are ejected, the cluster core shrinks again on a 
time scale of ^ time units. 

About one-tenth of the known globular clusters in the 
Galaxy have vanishingly small cores and are inferred to be in 
a state of core collapse (Harris 199^. We note here the large 
number (^ 10) of GCs in Figure[2with small or zero core 
radii; this may indicate that a much larger fraction 80%) 
of the LMC globular clusters are on their way to core col- 
lapse. We predict that these clusters have lost most or all of 
their black holes, while the ~ two old clusters in Figure^with 
substantial cores still contain a few stellar mass black holes in 
their cores. We find no indication from their structural pa- 
rameters that these two clusters differ systematically from the 
other clusters in Figure 1 . 

Figure|3lshows the evolution of the core radii in these simu- 
lations. Computation of was based on its standard definition 
as the projected radius at which the surface density falls to 1 /2 
of its central value. Projected densities were computed via a 
kernel estimator ( Merritt & Tremblav 1994: Merritt2004). To 
reduce the noise, values of from all experiments with the 
same 7 and with Nbh =10 were averaged together Figure 3 
shows that core sizes increase roughly as the logarithm of the 
time, consistent with the time dependence of the upper enve- 
lope of Figure 1 , and reach values at the end of the simulations 
of ^ 10% of the half-mass radius. 

Based on equation (6) and the discussion following, the 
conversion factor from model time units to physical time units 
is approximately 



8.9 X \{fyrM^^'^a^QN^,.s 



(7) 
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with A^^ 5 the number of stars in the GC in units of 10^. This 
scaHng was used to plot three curves in Figure 1 : with M5 = 
N^^i = 2, flio = 0.5 (bottom), Mg = A^*,5 = 0.5, flio = 1 (middle), 
andMg = A^*.5 = 1, oio = 2.5 (top). The curves in Figure 1 were 
taken from the experiments with 7=1; the experiments with 
7 = 1.5 and 2 give similar results (note that R^/a varies by a 
factor ~ 2 from 7 = 1 to 7 = 2, hence rc/Re varies less than 
rc/a in Figure 3). The logarithmic time dependence of the 
upper envelope of the distribution is well reproduced, and 
with appropriate (and reasonable) scaling, points below the 
envelope can also be matched. As noted above, the smaller 
core radii that begin to appear in SMC/LMC clusters with t > 
10^ yr are plausibly due to evolution toward core collapse in 
these clusters, as seen also in some of the simulations. 
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Fig. 3. — Evolution of the core radius, defined as the radius at which 
the projected density falls to one-half of its central value. Each curve is 
the average of the various experiments at the specified value of 7, with 
7 = 0.5 (top), 7=1 (middle) and 7 = 2 (bottom). Vertical axis is in units 
such that the Dehnen-model scale length a = 1; see text for conversion 
factors from a to 

4. DISCUSSION 

The core formation mechanism proposed here could begin 
to act even before the most massive stars had evolved into 
black holes. Evolution times for 2OM0 stars are ~ 8 Myr 
(^challeret al. 1992). The earliest phases of core formation, 
r < 10^ yr, would therefore be driven by the accumulation of 
massive stars rather than by their remnants. Figure 1 shows 
possible evidence of core growth on time scales < 20 Myr 
in a few clusters. In this context it is interesting to mention 



the so-called "young dense star clusters." These clusters have 
ages ^ 10 Myr, sizes ~ 1 pc, and c ontain ^ 10^ stars . Well- 
known examples are NGC2070 (Brandl et al." 1996V NGC 
3603 (Vrba 2000), and Westerlund 1 (Brandl etal. 1999). 
All of these young clusters have small but distinct cores. 
The young cluster R136 in 30 Doradus (r ~ 5 Myr) shows 
clear evidence of m ass segregation among the brightest stars 
jBrandl etalJl996h . 

The mechanism described here is similar to core formation 
by a binary supermassive black hole in a galactic nucleus via 
the gravitational slingshot (Ouinlan 1996). The latter pro- 
cess produces cores with masses ~ a few times the binary 
mass, assuming that the binary separation decays all the way 
to the point that coalesc ence by gravitational wave emission 
can ensue (Merritt 2003). If the decay stalls at a larger separa- 
tion, the displaced mass will be smaller. It is currently uncer- 
tain h ow often the decay would stall (Milosavlievic & MerritJ 
l2003h . Core formation by a population of massive remnants 
also displaces a mass that is a few times the total mass in black 
holes (Figure 2), and because the smaller black holes are freer 
to move about, there is less prospect of stalling due to a lo- 
cal depletion of stars. In galactic nuclei, the imprints left on 
the stellar distribution by the clustering of stellar-mass black 
holes were probably long ago erased by the growth of the su- 
permassive black hole, by the formation and decay of binary 
supermassive black holes during galaxy mergers, and by star 
formation. 

A speculative application of these results is to cores 
formed at the centers of dark matter halos by the clus- 
tering of Population III remnants in the early universe 
((Volonteri, Haardt, & Madau 2003). The latter are believed 
to contain at least one- half the mass of their stellar pro geni- 
tors when m > 25OM0 d Frver. Wooslev. & He gedl2001 ). and 
the cosmological density of remnants may be similar to that of 
the supermassive black holes presently observed at the centers 
of galaxies (Madau & Rees 2001). It follows that the Popu- 
lation III remnants could create cores of appreciable size, ;/ 
a number of them can accumulate in a single halo at a given 
time, and if the time for their orbits to decay is shorter than 
the time between halo mergers. Both propositions will require 
further investigation. 
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